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We present the numerical Quantum Monte Carlo results for the ground state energy of circular 
quantum dots in which Rashba spin-orbit iteraction is present. Diffusion Monte Carlo with spin 
propagation is applied in order to treat the spin-orbit interaction correctly, following previous work 
done in the field of the two-dimensional electron gas. Together with ground state energies, also 
numerical results for density and spin-density profiles are given. 
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I. INTRODUCTION 

Quantum dots have become in the last years very im- 
portant in several fields, from spintronics to solar pan- 
els and even biological applications. They represent also 
possible candidates for qbits realization in view of quan- 
tum computation. The reason for such a success is that 
they are relatively easy to produce and their discrete en- 
ergy levels^ can be engineered depending on specific 
needs acting over the intensity and shape of the confin- 
ing potential. Though a quantum dot is indeed a quasi 
zero-dimensional structure, in which electrons are con- 
fined along all three dimensions, several possible real- 
izations of quantum dots exist. In this work we con- 
sider in particular quantum dots obtained as finite size 
quasi two-dimensional structures at semicondutor inter- 
faces. A semiconductor heterostructure provides a very 
narrow (001) quantum well from which the quasi two- 
dimensional character raises. The lateral confinement 
may then be obtained applying an electrostatic poten- 
tial or using etching tecniques. Modern fabrication tech- 
niques allow also the control of the number of the con- 
fined electrons, which can be very low. 

Due to a structure inversion asymmetry of the ma- 
terial an electric field generates inside of the quantum 
well along the z direction causing the appearance of the 
Rashba interaction^: 
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Here Pi is the momentum of the i — th electron, and 
of and of are the Pauli matrices acting over the spin 
of particle i. This spin-orbit interaction is particularly 
interesting due to its tunability^r— : acting on the gate 
voltage it is possible to control the coupling constant 
A. This in principle gives the possibility of modifying 
the system spin and the energy levels through a easy 
to operate electric potential. Spin-orbit effects have a 
particular interest in confined systems, due to the pos- 
sible appearance of features like non obvious spin densi- 
ties and for their effects on energy levels and shell struc- 
tures. Several Monte Carlo simulations have already 



been done in different conditions^— , proving the valid- 
ity and the accuracy of the method. Other techniques 
have also been applied, such as Hartree-Fock, Den- 
sity Functional-based methods^ like LDA and LSDA. 
Also exact diagonalizatio n 1 b 12 was used, but mainly for 
very small numbers of electrons. Quantum dots with 
Rashba interaction were studied in particular in absence 
of Coulomb potentiate and recently a DFT calculation 
has been done also taking the Coulomb interaction into 
account^. What we propose here is an ab-initio Monte 
Carlo study, in which the spin-orbit interaction is taken 
into account by a spin-operator dependent propagator. 



II. METHOD 

A. Hamiltonian 

The system studied in this paper is a circular quantum 
dot with different numbers of electrons trapped into a 
static parabolic electric potential of the form V con f (r) = 

The Hamiltonian for 



-/2, where 



such a system with N electrons in presence of the Rashba 
spin-orbit interaction can be written as 
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Here e and m* respectively are the effective dielectric 
constant and electron mass in the semiconductor. About 
the confining potential strenght u and the Rashba in- 
teraction coupling constant A different values have been 
used in our calculations in order to obtain a more com- 
plete view of the contributions of the two interactions to 
the system ground state. In the following we will make 
use of effective atomic units in order to simplify notation, 



defining h = e /e 



1. The lenght unit therefore 



is the effective Bohr radius Oq = aoe/m* and the energy 
unit the effective Hartree H* = H ■m* / (m e e 2 ). For GaAs 
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quantum dots, we take e = 12.4 and to* = 0.067to 6 , yeld- 
ing H* = 11.86meV and a* = 97.93A This Hamilto- 
nian, due to the spin-orbit term, does not commute with 
the angular momentum operator, but can be proven to 
commute with the z component of the total momentum 

Jz- 



B. Diffusion Monte Carlo 

Diffusion Monte Carlo is a very accurate method used 
for the study of the ground state properties of many 
body systems. It is based on projection in imaginary 
time of an initial wavefunction, exploiting the fact that 
higher energy components of the initial wavefunction will 
tend to be exponentially suppressed respect to the ground 
state. In order to obtain imaginary time evolution of a 
wavefunction an imaginary time propagator is needed. 
Though an exact form for the propagator is unknown, 
Trotter's formula is employed in order to write the prop- 
agator as a product of exponentials of the different terms 
of the Hamiltonian. This is correct only to second or- 
der in the propagation time. The problem is overcome 
by applying a short time propagation repeatedly in order 
to obtain a long enough projection. The DMC method 
samples the ground state wavefunction (multiplied by a 
trial wavefunction when importance sampling is used) 
with walkers, i.e. points in the phase space moving dur- 
ing propagation. The kinetic term of the Hamiltonian 
produces a gaussian propagator, which is used in DMC 
as a sampling probability for walkers displacements. The 
propagator factor coming from the potential is instead 
usually seen simply as a weight for the walkers. In case 
the potential contains terms which only depend on space 
coordinates this is in principle not a problem. 
For momentum depending potentials anyway, this ap- 
proach cannot be applied. Being walkers points in the 
coordinate space they don't have a definite momentum, 
and in almost all interesting cases also the trial wave- 
function is not a momentum eigenstate. The Rashba po- 
tential, being due to a spin-orbit-like interaction, shows 
indeed dependence on the particles momenta and there- 
fore needs to be treated in a particular way. The basic 
idea (for a detailed description of the method seei£) is 
dividing the propagator into three factors, containing re- 
spectively the kinetic energy, the spin-orbit potential and 
the remaining parts of the Hamiltonian (which then only 
depend on space coordinates). One can choose to apply 
the factor containg the Rashba potential right after the 
gaussian free particle propagator. With this procedure 
the derivative terms inside of the spin-orbit interaction 
are turned into terms with space coordinates and spin 
operators dependence, having the meaning of a spin ro- 
tation. The third factor can finally be applied in the 
usual way, i.e. as a weight. The full propagator for the 



quantum dot in this case can be written as: 

G(R, R', At) = e -[ WRR W*)-*o-^]Ar x (3) 
e -^EL(A^-A<^)Ar Gfl(RiR ' !AT) . 

With R and R' we respectively indicate the set of new 
and old space coordinates of the N electrons, with Vcoui 
the Coulomb potential and with Arj (j = x, y) the x and 
y components of the difference between the new and the 
old coordinates of the i-th electron. Go indicates the free 
particle propagator and At the imaginary time step for 
the evolution. Notice that the propagator has already 
been renormalized by the factor cxp(£'oAr), where Eq is 
the ground state energy, in order to let the ground state 
weight remain finite. 

C. Wavefunction 

Though the quantum dot with parabolic confinement 
and Rashba interaction can be described in absence of 
Coulomb potential by a sum of single particle Hamilton- 
inans, its analytical solution has nevertheless not been 
found yet. The problem one encounters when trying to 
diagonalize such Hamiltonians over a basis of harmonic 
oscillators, is an infinite set of coupled equations. What 
has been so far proposed is an analytical very accurate 
approximation and a numerical diagonalization over two 
dimensional harmonic oscillator states^. Due to the fact 
that [H, J z ] — Q one expects to get better results for 
the ground state energy using a trial wavefunction which 
is a J z eigenstate since such a state would have the re- 
quired symmetry. Because the single particle orbitals 
deriving from the two methods mentioned above mix dif- 
ferent harmonic oscillator states, it is difficult building 
with them a J z eigenstate. On the other hand one can 
easily provide J z eigenvectors with simple combinations 
of two dimensional harmonic oscillator orbitals but such 
wavefunctions will be both eigenstates of S z and L z and 
therefore will be eigenstates of operators which do not 
commute with the Hamiltonian. From our Monte Carlo 
calculations wc observed that the wavefunctions obtained 
from simple combinations of harmonic oscillators always 
give the lowest results for the ground state energy, there- 
fore we adopted them obeying to the variational princi- 
ple. 

Because of the spin rotating factor contained in the imag- 
inary time propagator just shown in the previous section, 
walkers will not only change their space coordinates dur- 
ing propagation, but also their spins. This implies that 
it is not possible separating spin-up and spin-down elec- 
trons in two Slater determinants: only one determinant 
must be used containing all electrons. Besides this it was 
necessary dealing with complex numbers because though 
single particle orbitals could be real, spin components in 
general are not. For these reasons the trial wavefunction 
is complex and its phase is not constant. This means 
that using the fixed node approximation is not possible. 
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uj = 0.28 



N e i L S A = A = 0.1 A = 0.35 A = 0.7 

2 0.51081(3) 0.5118(3) 0.5120(3) 0.5153(3) 

3 1 1/2 0.7446(2) 0.746(1) 0.748(1) 0.753(1) 

4 0.94520(5) 0.944(1) 0.946(1) 0.951(1) 

4 1 0.92863(5) 0.934(1) 0.934(1) 0.939(1) 

5 1 1/2 1.1067(2) 1.117(1) 1.118(1) 1.127(1) 

6 1.2666(4) 1.282(1) 1.285(1) 1.289(1) 

7 2 1/2 1.4337(4) 1.457(1) 1.463(1) 1.493(3) 

ui = 1.78 



N el L S 



A = A = 0.1 A = 0.35 A = 0.7 



2 2.4820(6) 2.4805(6) 2.4803(5) 2.4815(5) 

3 1 1/2 3.469(1) 3.4663(4) 3.4669(5) 3.485(5) 

4 4.200(2) 4.1968(4) 4.1969(4) 4.1975(4) 

4 1 4.150(1) 4.151(1) 4.148(1) 4.150(1) 

5 1 1/2 4.784(1) 4.7832(5) 4.7839(5) 4.786(1) 

6 5.332(2) 5.326(1) 5.326(1) 5.329(1) 

7 2 1/2 6.003(1) 6.004(1) 6.003(1) 6.012(2) 

TABLE I: Diffusion Monte Carlo results for the ground state 
energy per particle for different numbers of electrons and val- 
ues of the Rashba coupling constant A. Results are reported 
for two different confinement potential streghts. Energies are 
expressed in effective Hartrees (H*). 



where n(N) is the chemical potential of a N electrons 
quantum dot. In both cases it is possible noticing that 
the behavior of E a dd as a function of the electrons num- 
ber becomes consistently steeper with increasing A. At 
lu = 1.78 comparison is possible with data from LSDA 
calculations^. Such comparison gives resonable agree- 
ment for A = but show different dependence on the 
Rashba coupling constant. Deviations become bigger 
increasing A since LSDA predicts a change in the shell 
structure. 
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Our choice has been the Fixed Phase approximatior 
which is based on the idea of finding the lowest energy 
state having the same phase as the trial wavefunction. 



FIG. 1: (color online) Addition energy for u = 0.28. A = 
are represented by the line with dots (black) ,A = 0.1 dashes 
(red) ,A = 0.35 dot-dash (green), A = 0.7 double dot-dash 
(blue). Data are given in effective atomic units. 



III. RESULTS 

The results contained in this paper concern the ground 
state of the system for variable numbers of electrons. The 
energy was calculated for two values of the parabolic con- 
fining potential and for three values of the Rashba inter- 
action straight in order to give an extensive and sys- 
tematic study of the system. Studying the system under 
different conditions is meant to partially account for the 
tunability of the two potentials. 

The results reported in table Hill show the different be- 
havior of the system depending on the value of uj. For 
smaller confinement it is clear how the spin orbit interac- 
tion always increases the ground state energy. Moreover 
the energy always shows the same increasing trend with 
increasing A. Something different happens for higher val- 
ues of the parabolic confinement. Though, due to error- 
bars it is not possible giving a detailed picture, ground 
state energies can be lower in presence of the spin orbit 
interaction, in particular for low values of the Rahsba 
coupling constant. They tend in any case to increase for 
high Rashba potential strenghts. 

In Fig. 1 1121 the addition energies for the system are re- 
ported respectively for w = 0.28 and u> = 1.78. The 
addition energy is defined as 



E, 



dd 



fi(N + l)-fx(N) 
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FIG. 2: (color online) Addition energies for dots with con- 
finement strenght u = 1.78 A = dots (black), A = 0.1 dashes 
(red), A = 0.35 dot-dash (green), A = 0.7 double dot-dash 
(blue). Data are given in effective atomic units. 

Quantum dots are known to show interesting density 
profiles, changing shape with the electrons number—. In 
particular, a 6 electrons dot shows a little hole in the 
density profile in its center. Such distributions are due 
to the interplay between the Coulombic repulsion and 
the kinetic energy within the confining potential. From 
Fig. |3] it is possible noticing how the Rashba interaction 
is effective in modifying the density of the dot, render- 
ing the central hole shallower. As one might expect the 
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spin orbit interaction does not allow the occupation of S z 
eigenstates since [H, S z ] is non zero once the interaction 
is switched on. This implies a different effect on average 
of the Coulomb potential and at the same time a modi- 
fication of the kinetic energy. 

Within the QMC approach it is possible sampling, apart 
from the density, also the spin density of the system, also 
studied within the DFT approach^. This correspond to 
the following expectation value: 



Pj( r ) 



< 



> 



<0(r')|0(r') > 



(5) 



with j = x,y,z. From Monte Carlo simulations such 
quantities in the two electrons case always gave fluctu- 
ating results, though always showing the same kind of 
structure in all cases (x,y,z). Such structures appeared 
with different orientations and peaks heights. This sug- 
gests that such spin densities are more related to corre- 
lation effects and should be connected to some two body 
operator. The results reported in Fig. [5] show the av- 
erage of 150 spin density results rotated in order to ob- 
tain alignment of the common structure along the x axis. 
From such figure it is possible noticing how the two dif- 
ferent spin components tend to arrange in opposite sides 
of the dot. The structures of Fig. [4] are common to all 
three spin density components {x,y,z) which give the 
same average results within errorbars (not reported in 
the graph, but of the order of 20 percent). 




FIG. 4: Rotation-averaged Y Spin Density of a 2 electons dot 
with confinement strenght u — 1.78 in presence of Rashba 
interaction A = 0.1. Data are reported in effective atomic 
units (do for lenghts) and (ao)~ 2 for the density. 




FIG. 5: The plot shows the X Spin Density of a 3 electons 
dot with confinement strenght ui = 0.28 in presence of Rashba 
interaction A = 0.1. Data are reported in effective atomic 
units (dp for lenghts) and (aj5)~ 2 for the density. 
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FIG. 3: (color online) Density of a 6 electons dot with confine- 
ment strenght to — 0.28 in presence of and without Rashba 
interaction. Triangles (blue) refer to results with A = while 
circles (green) represent results at A = 0.35. Data are given 
in effective atomic units. 

Calculations for the spin-density were done also for a 3 
electrons quantum dot with lo — 0.28 and A = 0.1. While 
in the 2 electrons system a non rotated average would 
give results consistent with zero spin-density within the 
errorbars, the 3 electron system shows interesting spin 
patterns and non zero magnetization. This is due to 
the fact that, while a two electron dot has a closed shell 
configuration, the addition of a third electron causes the 
partial occupation of the second shell, corresponding to 



non zero angular momentum and spin. From fig. I5I6I7I 
it is clear how the z spin density has a structure which 
strongly differs from the x and y spin densities. In par- 
ticular, within errorbars, which are of the order of 10 
percent, the z spin density has circular symmetry, while 
the x has a structure very similar to the y spin density, 
though having orthogonal orientation. It is interesting to 
mention how such structures show a general resemblance 
to the analytical results one would obtain in absence of 
Coulomb interaction using approximate solutions i 19 ' 20 . 
Since the trial wave function employed in QMC calcula- 
tions was a combination of gaussian orbitals and carried 
no information about the effects of the spin orbit inter- 
action on the system, we stress that such peculiar struc- 
tures are only due to the spin orbit propation given by 
the Green's function |QJ. 
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IV. CONCLUSIONS 




FIG. 6: Y Spin Density of a 3 electons dot with confinement 
strenght to — 0.28 in presence of Rashba interaction A = 0.1. 
Data are reported in effective atomic units (a£ for lenghts) 
and (do) -2 for the density. 



This study represents an application of the DMC 
method with spin-orbit propagator— in the direction of 
confined electrons system. It opens therefore the way to 
further applications of Diffusion Monte Carlo to nanosys- 
tems and molecules in presence of spin-orbit interactions 
which should not necessarily coincide with the Rashba 
interaction. The different structure of the energy levels 
depending on the confining potential suggests the possi- 
bility of having more freedom in the engineering of energy 
levels. Though, due to errorbars a full characterization of 
the energy levels structure of the system is not yet possi- 
ble, the results presented show nevcthclcss an interesting 
and non trivial dependence on the two parameters A and 
lu. The method allowed furthermore the study of spin 
orbit effects on the electron density and spin density. 
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